/*
 * Copyright (C) 2009-2010 Andre Schulz, Florian Jung, Sebastian Hartte,
 *						   Daniel Trick, Christan Wojek, Konrad Schindler,
 *						   Jens Ackermann, Michael Goesele
 * Copyright (C) 2008-2009 Christopher Evans <chris.evans@irisys.co.uk>, MSc University of Bristol
 *
 * This file is part of SURFGPU.
 *
 * SURFGPU is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 8
 * SURFGPU is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with SURFGPU.  If not, see <http://www.gnu.org/licenses/>.
 */

#include <cstdio>
#include <cstdlib>

#include <cuda_runtime.h>
#include <cudpp/cudpp.h>

#include "integral.h"
#include "convertRGB2GrayFloatGPU.h"
#include "transposeGPU.h"

//! Computes the integral image of image img.  Assumes source image to be a
//! 32-bit floating point.  Returns cudaImage of 32-bit float form.
cudaImage* IntegralMat(cv::Mat *src)
{
	unsigned int *d_rgb_img;
	float *d_gray_img;
	float *d_int_img;
	float *d_int_img_tr;
	float *d_int_img_tr2;

	// Allocate device memory
	int img_width = src->cols;
	int img_height =  src->rows;
	size_t rgb_img_pitch, gray_img_pitch, int_img_pitch, int_img_tr_pitch;
	
#ifdef DEBUG  
  printf("width: %d height: %d ",img_width,img_height);
  fflush(stdout);
#endif    
  cutilSafeCall( cudaMallocPitch((void**)&d_rgb_img, &rgb_img_pitch, img_width * sizeof(unsigned int), img_height) );
	cutilSafeCall( cudaMallocPitch((void**)&d_gray_img, &gray_img_pitch, img_width * sizeof(float), img_height) );
	cutilSafeCall( cudaMallocPitch((void**)&d_int_img, &int_img_pitch, img_width * sizeof(float), img_height) );
	cutilSafeCall( cudaMallocPitch((void**)&d_int_img_tr, &int_img_tr_pitch, img_height * sizeof(float), img_width) );
	cutilSafeCall( cudaMallocPitch((void**)&d_int_img_tr2, &int_img_tr_pitch, img_height * sizeof(float), img_width) );

	// Upload RGBA image to GPU
	cutilSafeCall( cudaMemcpy2D( d_rgb_img, rgb_img_pitch,
								 src->data, img_width * sizeof(unsigned int),
								 img_width, img_height,
								 cudaMemcpyHostToDevice ));
	
  // Convert RGB to gray
	convertRGB2GrayFloatGPU(
		d_gray_img, gray_img_pitch,
		d_rgb_img, rgb_img_pitch,
		0x00ff0000, 0x0000ff00, 0x000000ff,
		16, 8, 0,
		img_width, img_height);

	// Setup cudpp multiscan plans for computing the integral image
	CUDPPHandle mscan_plan, mscan_tr_plan;

	CUDPPConfiguration cudpp_conf;
	cudpp_conf.op = CUDPP_ADD;
	cudpp_conf.datatype = CUDPP_FLOAT;
	cudpp_conf.algorithm = CUDPP_SCAN;
	cudpp_conf.options = CUDPP_OPTION_FORWARD | CUDPP_OPTION_INCLUSIVE;

	CUDPPResult res = cudppPlan(&mscan_plan, cudpp_conf,
		img_width, img_height, gray_img_pitch / sizeof(float));
	if (res != CUDPP_SUCCESS)
	{
		std::printf("ERROR: Couldn't create CUDPP plan!\n");
		std::exit(EXIT_FAILURE);
	}

	res = cudppPlan(&mscan_tr_plan, cudpp_conf,
		img_height, img_width, int_img_tr_pitch / sizeof(float));
	if (res != CUDPP_SUCCESS)
	{
		std::printf("ERROR: Couldn't create CUDPP plan!\n");
		std::exit(EXIT_FAILURE);
	}

	cudppMultiScan(mscan_plan, d_int_img, d_gray_img, img_width, img_height);
	transposeGPU(d_int_img_tr, int_img_tr_pitch,
		d_int_img, int_img_pitch,
		img_width, img_height);
	cudppMultiScan(mscan_tr_plan, d_int_img_tr2, d_int_img_tr, img_height, img_width);
	transposeGPU(d_int_img, int_img_pitch,
		d_int_img_tr2, int_img_tr_pitch,
		img_height, img_width);

	cudppDestroyPlan(mscan_plan);
	cudppDestroyPlan(mscan_tr_plan);
	cutilSafeCall( cudaFree(d_rgb_img) );
	cutilSafeCall( cudaFree(d_gray_img) );
	cutilSafeCall( cudaFree(d_int_img_tr) );
	cutilSafeCall( cudaFree(d_int_img_tr2) );

	cudaImage *img = new cudaImage;
	img->width = img_width;
	img->height = img_height;
	img->data = reinterpret_cast<char*>(d_int_img);
	img->widthStep = int_img_pitch;

	return img;
}


cudaImage* Integral(IplImage *src)
{
	unsigned int *d_rgb_img;
	float *d_gray_img;
	float *d_int_img;
	float *d_int_img_tr;
	float *d_int_img_tr2;

	// Allocate device memory
	int img_width = src->width;
	int img_height = src->height;
	size_t rgb_img_pitch, gray_img_pitch, int_img_pitch, int_img_tr_pitch;
	cutilSafeCall( cudaMallocPitch((void**)&d_rgb_img, &rgb_img_pitch, img_width * sizeof(unsigned int), img_height) );
	cutilSafeCall( cudaMallocPitch((void**)&d_gray_img, &gray_img_pitch, img_width * sizeof(float), img_height) );
	cutilSafeCall( cudaMallocPitch((void**)&d_int_img, &int_img_pitch, img_width * sizeof(float), img_height) );
	cutilSafeCall( cudaMallocPitch((void**)&d_int_img_tr, &int_img_tr_pitch, img_height * sizeof(float), img_width) );
	cutilSafeCall( cudaMallocPitch((void**)&d_int_img_tr2, &int_img_tr_pitch, img_height * sizeof(float), img_width) );

	// Upload RGBA image to GPU
	cutilSafeCall( cudaMemcpy2D( d_rgb_img, rgb_img_pitch,
								 src->imageData, src->widthStep,
								 img_width * sizeof(unsigned int), img_height,
								 cudaMemcpyHostToDevice ));

	// Convert RGB to gray
	convertRGB2GrayFloatGPU(
		d_gray_img, gray_img_pitch,
		d_rgb_img, rgb_img_pitch,
		0x00ff0000, 0x0000ff00, 0x000000ff,
		16, 8, 0,
		img_width, img_height);

	// Setup cudpp multiscan plans for computing the integral image
	CUDPPHandle mscan_plan, mscan_tr_plan;

	CUDPPConfiguration cudpp_conf;
	cudpp_conf.op = CUDPP_ADD;
	cudpp_conf.datatype = CUDPP_FLOAT;
	cudpp_conf.algorithm = CUDPP_SCAN;
	cudpp_conf.options = CUDPP_OPTION_FORWARD | CUDPP_OPTION_INCLUSIVE;

	CUDPPResult res = cudppPlan(&mscan_plan, cudpp_conf,
		img_width, img_height, gray_img_pitch / sizeof(float));
	if (res != CUDPP_SUCCESS)
	{
		std::printf("ERROR: Couldn't create CUDPP plan!\n");
		std::exit(EXIT_FAILURE);
	}

	res = cudppPlan(&mscan_tr_plan, cudpp_conf,
		img_height, img_width, int_img_tr_pitch / sizeof(float));
	if (res != CUDPP_SUCCESS)
	{
		std::printf("ERROR: Couldn't create CUDPP plan!\n");
		std::exit(EXIT_FAILURE);
	}

	cudppMultiScan(mscan_plan, d_int_img, d_gray_img, img_width, img_height);
	transposeGPU(d_int_img_tr, int_img_tr_pitch,
		d_int_img, int_img_pitch,
		img_width, img_height);
	cudppMultiScan(mscan_tr_plan, d_int_img_tr2, d_int_img_tr, img_height, img_width);
	transposeGPU(d_int_img, int_img_pitch,
		d_int_img_tr2, int_img_tr_pitch,
		img_height, img_width);

	cudppDestroyPlan(mscan_plan);
	cudppDestroyPlan(mscan_tr_plan);
	cutilSafeCall( cudaFree(d_rgb_img) );
	cutilSafeCall( cudaFree(d_gray_img) );
	cutilSafeCall( cudaFree(d_int_img_tr) );
	cutilSafeCall( cudaFree(d_int_img_tr2) );

	cudaImage *img = new cudaImage;
	img->width = img_width;
	img->height = img_height;
	img->data = reinterpret_cast<char*>(d_int_img);
	img->widthStep = int_img_pitch;

	return img;
}

